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We consider different Markovian embedding schemes of non-Markovian stochastic processes that 
are described by generalized Langevin equations (GLE) and obey thermal detailed balance under 
equilibrium conditions. At thermal equilibrium superdiffusive behavior can emerge if the total 
integral of the memory kernel vanishes. Such a situation of vanishing static friction is caused by 
a super-Ohmic thermal bath. One of the simplest models of ballistic superdiffusion is determined 
by a bi-exponential memory kernel that was proposed by Bao [J.-D. Bao, J. Stat. Phys. 114, 
503 (2004)]. We show that this non-Markovian model has infinitely many different 4-dimensional 
Markovian embeddings. Implementing numerically the simplest one, we demonstrate that (i) the 
presence of a periodic potential with arbitrarily low barriers changes the asymptotic large time 
behavior from free ballistic superdiffusion into normal diffusion; (ii) an additional biasing force 
renders the asymptotic dynamics superdiffusive again. The development of transients that display 
a qualitatively different behavior compared to the true large-time asymptotics presents a general 
feature of this non-Markovian dynamics. These transients though may be extremely long. As a 
consequence, they can be even mistaken as the true asymptotics. We find that such intermediate 
asymptotics exhibit a giant enhancement of superdiffusion in tilted washboard potentials and it 
is accompanied by a giant transient superballistic current growing proportional to t a ° ff with an 
exponent a e ff that can exceed the ballistic value of two. 

PACS numbers: 05.40.-a, 82.20.Uv, 87.16.Uv 



I. INTRODUCTION 

The subject of anomalous diffusion has become increas- 
ingly popular and important in the last years with a num- 
ber of papers growing faster than linearly in time with 
almost 500 papers published last year. There also is a 
large number of theoretical models that lead to anoma- 
lous diffusion such as as continuous time-random walks 
[1-4], including Levy flights and Levy walks [3, 5], re- 
lated fractional Fokker-Planck equations [5, 6] and (ordi- 
nary) Langevin equations in random subordinated time 
[7, 8], as well as (ordinary) Langevin equations with ad- 
ditive non- Gaussian Levy white noises [9, 10]. Moreover 
nonlinear Brownian motion with multiplicative Gaussian 
white noise [11, 12], as well as linear Boltzmann equation 
with scattering events being distributed in time according 
to a power law distribution [13, 14] may display anoma- 
lous diffusion. This list by far is not complete. Yet the 
quest for minimal and fundamental physical models has 
become ever more important. One of the fundamental 
approaches to anomalous diffusion [15-18] is provided by 
the generalized Langevin equation (GLE) [19-22] with a 
frictional memory kernel "f(t), reading: 



+ m / 7 (t - t')x(t')dt' + 



dV(x,t) 
dx 



= m , (i) 



where x(t) denotes the position of a particle of mass m. 
Here £(i) is a Gaussian zero-mean fluctuating force that 
at temperature T is related to the memory kernel by the 
fluctuation-dissipation relation [19] 



(C(t)((t / ))=k B Tm 1 (\t~t'\) 



(2) 



Remarkably, this model can be derived from a Hamil- 
tonian dynamics of a particle that bilinearly cou- 



ples with coupling constants Cj to a thermal bath 
of harmonic oscillators with masses rrii and frequen- 
cies U>i, H B ,mt(pi,qi,x) = {l/2)J2i{Pi/ m t + m i^"i{li ~ 

CiX I (rriiUif)] 2 } . The total effect of the bath oscillators, 
which are initially canonically distributed with H B - lnt at 
temperature T and fixed x = x(0), is characterized by 
the bath spectral density 



J(w) = 



2 t— 1 rrijUJ. 



■6(u - w»). 



(3) 



It is related to the power spectral density of the fluctu- 
ating force 



(C(r)C(0))e- 



'dr 



(4) 



via S(w) = 2k B TJ(uj)/uj [19, 22]. This in general leads to 
a non-Markovian process of the particle dynamics with 
linear memory friction and Gaussian fluctuating force. 
Moreover, due to the fluctuation-dissipation relation (2) 
it is compatible with thermal equilibrium in confining, 
time-independent potentials and encompasses a whole set 
of physically meaningful models characterized by differ- 
ent bath spectral densities J(uS). 

In the absence of any potential the variance of the par- 
ticle's position will grow with time. The law according 
to which the variance grows characterizes the nature of 
the resulting diffusion process as being sub-diffusive if 
the growth of the variance is slower than linear. This 
happens if the static friction 7 = J* °° j(t')dt' diverges. 
Normal diffusion corresponds to a linear growth. It oc- 
curs if 7 is finite. Finally, if 7 vanishes the variance 
grows faster than linear and one speaks of superdiffu- 
sion. The presence of a nonlinear time-dependent force 
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f{x,t) = —dV(x,t)/dx modifies this simple picture in 
a complicated way depending on further details of the 
memory kernel and also on temperature. The qualitative 
behavior of the variance of the position is determined by 
the variance of the velocity, (Av 2 (t)}. The mean square 
displacement of position spreads according to normal dif- 
fusion, if the integral of the velocity variance over all 
times is hnite. On the other hand if this integral is zero, 
then the motion is anti-persistent and subdiffusive. If 
this integral diverges the spread of the position variance 
is superdiffusive. For free motion in the absence of a 
potential these criteria are equivalent to those for the 
memory kernel which in general fail in the presence of a 
potential. Since general analytical results are scarce and 
most likely nonexistent for nonlinear and time-dependent 
forcing, the reliability of numerical simulations has be- 
come a key issue. Numerically tractable models can be 
obtained by approximating the given memory kernel by a 
finite sum of exponential functions. The according non- 
Markovian particle dynamics can then be obtained as the 
projection of a high-dimensional Markovian process onto 
the phase space of the particle spanned by the particle's 
coordinate and momentum p = mx. The dimensional- 
ity of the Markovian process is D — N + 2, where N 
is the number of exponentials in the sum approximating 
the memory kernel. The key point is that the corre- 
sponding Markovian dynamics can be propagated locally 
in time for very long time intervals by means of very 
reliable algorithms with a well controlled numerical pre- 
cision. Moreover, this way of thinking allows one to iden- 
tify the simplest models for the superdiffusive GLEs with 
minimal embedding dimensions D = 3 and D = 4. The 
case D — 4 corresponds to approximating the memory 
kernel by a difference of two exponentials, 

-f(t) = 7i exp(-M) - 72 exp(-fc 2 *), (5) 

such that 71/fci =72/^2 implying to vanishing static fric- 
tion 7 and 71 > 72. The latter condition amounts to the 
fact that the memory kernel is proportional to the auto- 
correlation function of the fluctuating force ( and hence 
must be non-negative semidcfinite. This bi-cxponential 
model was proposed by Bao [23]. It corresponds to a 
super-Ohmic spectral density of thermal bath oscillators, 
J(uj) cx uj 3 at low frequencies and describes the coupling 
of a particle to three-dimensional lattice phonons. It 
therefore models the diffusion of an impurity in a crys- 
tal. We will demonstrate that this model can be em- 
bedded in infinitely many ways. In the following we will 
study one of the simplest embeddings which is different 
from those used by Bao [23-25]. In the absence of any 
force f(x,t) the spreading of the particle's position dis- 
tribution is ballistic, (Ax 2 (t)) = D 2 t 2 , and hence super- 
diffusive. Here and in the following, the expectation (...) 
refers to an ensemble average with respect to the fluc- 
tuating force £(i) and an initial distribution of position 
and momentum, x(0) and p{0), respectively. The velocity 
process is non-ergodic [26]. As a consequence, the bal- 
listic superdiffusion coefficient D 2 turns out to depend 



on the initial velocity distribution. However, as we shall 
see below this non-ergodic feature disappears in periodic 
potentials. Moreover, our numerics reveals that the bal- 
listic diffusion in tilted periodic potentials does neither 
depend on the initial velocity distribution, nor on the ini- 
tially non-equilibrium noise preparation. For this reason, 
we presume that the velocity process of the particle then 
becomes ergodic. 

The minimal, three-dimensional Markovian embedding 
of the GLE superdiffusion is achieved in the limit 71 — > 
00, ki — > 00, so that 70 = 71/fci = const = 72/^2- In 
this limit, the first exponential becomes a delta- function 
270 8(t). This case will be studied elsewhere. 

Unfortunately a non-Markovian Fokker-Planck-type 
equation (NMFPE) that corresponds to a GLE with a 
general, nontrivial potential is not known in spite of 
many years of search. The only exceptions are provided 
by (strictly) linear and parabolic potentials, where the 
corresponding NMFPEs were derived by Adelman [27] 
and Hanggi [28-30] for stable non-Markovian Brownian 
motion and, as well, for unstable non-Markovian dy- 
namics [31]; i.e. the Kramers problem of escape over 
a parabolic barrier [22]. Using the fact that [x(t),x — 
v(t)] is a two-component Gaussian process, which is ob- 
tained by a linear integral transformation of the Gaus- 
sian noise process £(t) the resulting NMFPE assumes 
the form of a time-dependent FPE. This FPE-structure 
of time-evolution for the single-time event probability of 
the non-Markovian process should then not be mistaken 
as an effective Markovian dynamics [28-30]. Notwith- 
standing these exceptional known cases cases of non- 
Markovian Gaussian dynamics, this lack of a generally 
closed NMFPE for nonlinear forces lends even more im- 
portance to the Markovian embedding approach. 

This paper is structured as follows. In Sec. II and 
Appendix I we detail the Markovian embedding proce- 
dure in a slightly more general way than has been used 
so far. The general results are illustrated with two differ- 
ent embeddings of one and the same superdiffusive GLE 
dynamics. In Sec. Ill, we present and discuss the results 
of stochastic simulations of superdiffusion under a con- 
stant bias and in a washboard potential for one of these 
embeddings. The issue of ergodicity of mean square dis- 
placement is discussed in Sec. IV. Conclusions arc drawn 
in Sec. V. 



II. METHOD 

The idea that we pursue here is to represent the non- 
Markovian stochastic dynamics of a single particle with 
(x,p) phase space as a projection of a multi-dimensional 
Markovian dynamics. It is well known that any GLE can 
be derived from the (Markovian) Hamiltonian dynamics 
of a particle coupled to a thermal bath of harmonic os- 
cillators. This Hamiltonian embedding though requires 
a large number of auxiliary degrees of freedom represent- 
ing the thermal bath. Here we look for an embedding 
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with a minimal number N of auxiliary variables which 
all together constitute a continuous Markovian process. 
A low embedding dimension is crucial for running nu- 
merical simulations which can become extensively time- 
consuming for a large N. 

We first rewrite the GLE (1) in terms of the phase 
space coordinates x and p = mx as 



x(t) =-p(t) 
m 



p(t) =f(x, t) - [ 7 (i - t')p(t')dt' + C(t) . 
Jo 



(6) 



The embedding involves a yet to be determined number 
N of auxiliary dynamical variables collected into a vector 
u(t) in terms of which the dynamics takes the following 
general form 

1 



x{t) =—p(t) 
m 

p(t) =f(x,t)+Tu(t) 

lt(t) = - p(t)f - At(t) + C£(i) 



(7) 



where ~g and ~r denote constant vectors of dimension N 
and A, and C are constant N x N matrices. The upper 
index T denotes the transpose of a vector or a matrix. 
Further, £ (t) is a vector of uncorrelated Gaussian white 
noises, 



(8) 



with N components. Integrating the equation for the 
auxiliary vector u(t) and substituting the result in the 
equation for the momentum p(t) one recovers the original 
GLE (6) only under special conditions, see Appendix A 
for details of the derivation. First, the memory kernel 
j(t) must satisfy 



7(*) 



— >T — At— > 



(9) 



Since the right hand side can in general be represented as 
a sum of N exponential functions exp(— Xit), i = 1 . . . N 
with the eigenvalues \ of the matrix A the embedding 
can be exact only if the memory kernel is of the same type 
[32] . But also other memory kernels such as algebraically 
decaying functions can be approximated by a finite sum 
of exponential functions, even with a relatively small ex- 
tra dimension N, and hence are amenable to Markovian 
embedding. 

Furthermore, the fluctuation-dissipation relation (2) im- 
poses restrictions on the matrices A, C, and the vectors 
~g and ~r. These restrictions are met if the embedding 
parameters satisfy the following two relations: 



G~g — rnksT^r 
CC T = AG + GA T 



(10) 
(11) 



which defines the constant N x N matrix G. 
However, for arbitrary initial values of the auxiliary vari- 
ables u(0), the fluctuation-dissipation relation will be 



obeyed only asymptotically. This means that the noise 
£(t) in Eq. (6) is initially nonstationary and becomes only 
gradually stationary in the course of time, see Appendix 
A, Eq. (A8). In order to guarantee the Gaussian nature of 
the random force the vector u(0) must also be Gaus- 
sian distributed, see Eq.(A3). Because the vector u(0) is 

independent of the vector of Gaussian white noises £ (t) 
and its first moment must vanish, it is sufficient to specify 
its covariance matrix 



<tT(0) ® 1I A (0)> = G . 



(12) 



It must coincide with G in Eq. (10), in order to have the 
fluctuation-dissipation relation (2) obeyed for all times, 
see Eq. (A8). 

Yet the conditions (9), (10) and (11) do not uniquely 
determine the enlarged process and actually leave room 
for an infinite variety of different processes leading upon 
reduction to the same generalized Langevin equation. 
Since some of the enlarged processes allow faster and 
more reliable numerical simulations than others there is 
a great interest in identifying computationally optimal 
embeddings. We further note that the relations (9), (10) 
and (11) are sufficient but not necessary conditions. The 
resulting embedding is more general than previous ones 
such as those proposed in Rcfs. [33, 34] which assume 
//• 

Before discussing a particular example we would like to 
emphasize that the stationarity of the fluctuating force 
and the fluctuation dissipation relation (2) are exactly 
implemented. 



A. Minimal model 

We now consider the simple class of models specified 
by the bi-exponcntial memory kernel (5). Under the con- 
dition of vanishing static friction, i.e. -fik 2 — 72&1, this 
memory kernel is specified by three independent param- 
eters that can be written as 

k 2 = 7i-72, v — k\ + k 2 and = k±k 2 ■ (13) 

Note that k is always a real parameter due to the positiv- 
ity constraint 71 > 72- In terms of this parameterization 
the Laplace transform 7(3) of the memory kernel becomes 



K 2 S 



7 (s)= / e - st 1 (t)dt= , 



(14) 



One can now easily calculate the power spectral density, 
see Eq. (4) , by connecting it to the friction kernel via the 
fluctuation-dissipation relation, see Eq. (2): 

S(uj) = 2k B Tm 7(i) cos(ujt)dt 
Jo 

= 2k B TmRe (j(iui)) (15) 
2k B TmK 2 vuj' 2 



(CO 2 — LO 2 ) 2 + V 2 UJ 
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It is interesting to note that the power spectral density 
has a maximum at the frequency cu a . 

There are several possibilities to realize Markovian em- 
bedding of this non-Markovian model. In the following 
we shall discuss two of them. 



1. First embedding 



A simple embedding is obtained by choosing: 



V CJ 

-uj a 



t T = («, 0) 



C = yj2mk B Tv ^ 
(ui(O)tij(O)) = mk B Tb l3 . 



1 




(16) 



(17) 



This choice leads to the following equations: 



x(t) =-p(t) 
m 

U\{t) = — Hp{t) — UU\{t) — UjQU 2 {t) 

u 2 (t) =u ui(t) , 



where £(t) is scalar Gaussian white noise. Moreover, this 
embedding also allows for complex parameters k\ — k 2 , 
providing the possibility to model oscillating real valued 
kernels 



7 (t) = K 2 e' ut/2 [cos{t^uj 2 -u 2 /A) 



(18) 



sin^u 2 - ^/4)] . 



This model corresponds to sharply peaked power spectral 
density S(ui), and bath spectral density J(u>). 



2. Second embedding 

An alternative way is to start with a diagonal matrix 
A. It would be tempting to also choose diagonal matri- 
ces C and G. This choice though always yields a linear 
combination of exponential functions with positive coef- 
ficients for the memory kernel (5) [35] and hence does not 
allow vanishing static friction. However, this goal can be 
achieved by means of the following choice of parameters 



involving non-diagonal matrices C and G 
A - [ fcl 

~ ' /,-, 



.91,2 = \ 7l,2 



fei + k 2 



ri,2 



C 



' k\ — k 2 

fci ^ h 
lh2 h + k 2 

\J 2mk B T 



(19) 



G 



where < c = 2\Jk\k 2 j (ki + k 2 ) < 1 is the correlation co- 
efficient of the covariance matrix G and k\ — k 2 (71/72) > 
k 2 . This choice is similar to the one in [23]. We note 
that the second embedding requires that the parameters 
k\ and k 2 must be real. Therefore it is not possible to 
model an oscillating kernel by this method. Our first em- 
bedding in Eq. (17) is, however, simpler and numerically 
more convenient since its numerical simulation requires 
less operations. For instance, only one stochastic vari- 
able has to be generated in the first embedding scheme, 
see Eq. (16), in contrary to two, needed in the second 
embedding scheme, see. Eq.(19). All this makes our first 
embedding scheme preferential. 



B. Dimensionless units 

For further studies, we transform Eq. (17) into dimen- 
sionless units by scaling momentum in terms of thermal 
momentum pt = ^/rnksT, expressing the distance in 
terms of a typical length scale some arbitrary length 
xo, which becomes the spatial period for periodic po- 
tentials (see below) and time in units of r = xq/vt- 
The auxiliary variables u\ and u 2 are scaled in units of 
uq = mxo/(KTg). The energy is scaled in units of k B T. 
This yields the equations of motion 



x =v 



V 



d ~ 

— V(x,t)+u 1 



u\ -- 
u 2 =Cj u\ 



k 2 v 



vu\ — Cjqu 2 + kv2z>£ 



(20) 



which were used in our simulations. Here, k = ktq, ujq = 
u; t , v = vtq. All results in the following figures are 
given in these dimensionless units. 



III. RESULTS 

The numerical results presented below were obtained 
using the standard stochastic Euler method [36]. A 
Mersenne Twister pseudo random number generator was 
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used to produce uniformly distributed random numbers 
which were transformed into Gaussian variables using 
Box-Muller algorithm [37]. Typically, an ensemble of 
n = 10 4 particles (or trajectories) was propagated in time 
with a fixed time step between At = 10~ 4 and 10 ~ 5 in 
most simulations to achieve (weak) convergence of en- 
semble averaged results. The use of double precision 
thus cannot be avoided and reliable numerics are very 
time consuming. All the particles were initially localized 
at x(0) = with the initial velocities sampled from some 
probability distribution. In most simulations we assumed 
this distribution to be sharply peaked at zero and as- 
cribed zero initial velocities to all the particles, although 
the thermal Maxwellian distribution was also used. The 
auxiliary variables Uj(0) were (mostly) sampled from the 
corresponding Gaussian distributions to achieve at the 
exact equivalence of the simulated Markovian dynamics 
to that of GLE, as described in the Section II. Method. 
Sometimes, we used also a different initial distribution 
of Ui(0) (all equal zero) in order to clarify the influence 
of the initially non-equilibrium noise preparation on the 
stochastic dynamics. In all cases, we denote the corre- 
sponding ensemble averages as (...) and specify the initial 
distributions if not obvious. 

Of central interest are the first moment (Ax(t)) and 
the variance (Ax 2 (t)} of the displacement 



Ax(t) = x(t) - x(0) = / dt'v(t'). (21) 
Jo 



Accordingly we have 



<As(t)> 



= f dt'(v(t' 
Jo 



)) 



(22) 



and 



(Ax 2 (t)) = ((Ax(t)-(Ax(t))) 2 ) 

d ft (23) 



=/7 

Jo Jo 



dtidt 2 C v (ti,t 2 ) 



where C v {t\,t 2 ) denotes the velocity fluctuation autocor- 
relation function 

C v (h,t 2 ) = ({v(h) (v(h)))(v(t 2 ) (v(t 2 )))) . (24) 

These quantities were estimated on the basis of averages 
over the ensemble of simulated particle trajectories. 

Of particular interest will turn out the question un- 
der which conditions the process of velocity fluctuations 
defined as the deviation of velocity from its mean value 
constitutes an ergodic process [38, 39]. The definition 
and main properties of an ergodic process are collected 
in Appendix B. 



A. Superdiffusion in presence of a constant bias 

First, we consider the Langevin dynamics (1) with an 
arbitrary memory kernel under a constant biasing force 



F, i.e with V(x, t) — —Fx. This special biased prob- 
lem is analytically solvable, cf. [15-17, 34], and there- 
fore provides a suitable test of our numerical simula- 
tions. The mean square displacement in this case does 
not depend on the external bias F and is given by: 
(Ax 2 (t)) = (x 2 (t)) - (x(t)) 2 becomes 

(Ax 2 {t)) = 2v% f H(t')dt' + [(v 2 (0)) - 4] H 2 (t),{25) 
Jo 

where we denote the thermal average of initial velocities 
by v\ = (v 2 (0)) T = k B T/m and 



H(t) = f K v (r)dT 
Jo 



(26) 



is the integral of the (normalized) equilibrium autocor- 
relation function of the velocity fluctuations which is de- 
fined as 



K v (t) = C v (t,0)/v%. . 



It has the Laplace-transform 



K v (s) 



s + 7(s) 



(27) 



(28) 



We note that the velocity fluctuations present a wide 
sense ergodic process if and only if the time average of 
K v (t) vanishes, i.e.: 



lim -H(t) = 0, 

t^oo t 



(29) 



see also Appendix B. For the mean displacement one 
obtains 



(Ax(t)) = — [ H{t')dt'. 
m Jo 



(30) 



If one chooses thermally distributed initial velocities, 
(w 2 (0)) = v%<, then the first and second moments 
of the displacement are connected by the fluctuation- 
dissipation theorem (FDT) 



(Ax(t)) 



F 



2k B T 



(Ax 2 (t)), 



(31) 



for any memory kernel. Notice that the mass of the parti- 
cle is not involved in Eq. (31). Provided that the velocity 
process is ergodic in the wide sense the second term on 
the right hand side of Eq. (25) can be neglected com- 
pared to the first term if time goes to infinity. Hence, for 
an ergodic velocity process the spreading of the particle 
position becomes independent of the initial velocity dis- 
tribution. In contrast, for a nonergodic process the sec- 
ond term can become comparable in magnitude or even 
dominant for large times. Then the influence of the ini- 
tial velocity distribution on the second moment of the 
position survives. This actually happens if the Laplace 
transform of the memory kernel, 7(s) approaches zero for 
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s — >■ proportionally to s or faster. In this case, the FDT 
(31) is not valid for (w 2 (0)) ^ v\, even asymptotically. 

As an example we consider the minimal model (5). Its 
Laplace transform indeed vanishes linearly with s — > 0, 
see Eq. (14). For the Laplace transform of the velocity 
correlation coefficient one obtains from Eq. (28) 



K v (s) 



w 2 



(32) 



Inverting the Laplace transform, one obtains for the ve- 
locity correlation coefficient 



K v (t) 



r/2 



(33) 



cosh(— — ^ t) 



^v 2 - 4w 2 - 4k 2 



sinh( 



sjv 2 - 4w 2 - 4k 2 



Note that with lim T 



K v (t) — luq/(w 2 + n 2 ) the equi- 



librium autocorrelation function of velocity fluctuations 
as well as its time average attain positive values. This 
confirms that the velocity process of the minimal model 
is nonergodic in the case of linear potentials. 

The mean square displacement of the position can be 
exactly evaluated by means of Eq. (25). We refrain from 
presenting the resulting lengthy expression and only com- 
pare the such derived exact result with the mean square 
displacement obtained from a simulation of the Marko- 
vian model via the first embedding for a particular set 
of parameters, see Fig. 1. The agreement between the 
analytical result and the simulation is very good. 

For short times the spreading of the mean square dis- 
placement of the position becomes 



(Ax 2 (t)) = (v\0))t 2 

+ k 2 [34-4(w 2 (0))] t 4 /12- 



0(t 5 ). (34) 



For a strictly vanishing initial velocity the contribution 
proportional to t 2 disappears and the diffusion initially 
becomes super-ballistic with (Ax 2 (t)) cx i 4 , see Fig. 1. 
Otherwise, the diffusion initially is ballistic. For large 
times ballistic diffusion results with (Ax 2 (t)) ~ DoX 2 . 
Due to the nonergodicity of the velocity process the bal- 
listic superdiffusion coefficient D2 depends on the initial 
distribution of velocities, 



D 2 = u 



T~7T 



! (0)> 



- 1 



(35) 



Fig. 1 also displays simulation results of the first 
embedding for vanishing initial values of the auxiliary 
variable, i.e. Wi(0) = u 2 (0) = 0, and also for ini- 
tial values from a Gaussian distribution with variance 



(«i(0) Ui (0)) 



< 
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numerical solution: (uj(0)u,(0)) = k 2 ^ 
numerical solution: ui(0) = 112(0) = * 2 ~ 
analytical solution 
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We recall that the latter choice 



FIG. 1: (Color online) The mean square displacement of 
the position (being independent of the bias F) as a function 
of time changes from a t 4 law at small times to a ballistic 
t 2 law. Comparison of the analytical solution according to 
Eq. (25) (solid line) and results from numerical simulations of 
the first embedding (square symbols) exhibit good agreement. 
A strongly deviating result is obtained if the auxiliary vari- 
ables of the first embedding initially assume vanishing values 
(dashed line). The other parameters k — 2, v — 3, ujo — 1 and 
v(0) — are the same for the displayed curves. The estimate 
of the mean square displacement is obtained by an average 
over an ensemble of 10 4 simulated trajectories. 



guarantees that the fluctuating forces are stationary and 
that they satisfy the fluctuation-dissipation relation (2). 
The mean square displacement resulting from zero initial 
auxiliary variables is remarkably different from that with 
the correct Gaussian distributed initial auxiliary vari- 
ables not only at short times but also for large times 
where it approaches normal instead of ballistic diffusion. 
The strong influence of the initial conditions even at large 
times is another consequence of the non-ergodicity of the 
velocity. In contrast, for an ergodic velocity process the 
long time behavior of the position mean square displace- 
ment has lost any memory on initial conditions. 



B. Superdiffusion in a washboard potential 

Next we consider the diffusion in a periodic washboard 
potential V(x, t) = —Vq cos(2-kx/x$) of spacial period xq. 
Here no analytical results are available, instead we per- 
formed numerical simulations of the first embedding. In 
Fig. 2 we compare the simulated mean square displace- 
ment as a function of time for different heights of the 
potential barriers separating neighboring periods of the 
potential. After a short initial period of fast growth the 
diffusion turns over in an intermediate ballistic behavior 
which eventually changes into normal diffusion. As far as 
one can say from the numerical simulations of finite dura- 
tion, normal diffusion always determines the asymptotic 
behavior. The onset time of normal diffusion though cru- 
cially depends on the magnitude of the potential barrier 
2Vq. The larger this barrier is the earlier normal diffusion 
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sets in. On the other hand for small barriers the ballistic 
regime extends over a large time before the asymptotic 
normal diffusion takes over. 




10" 1 10° 10 1 10 2 10 3 10 4 

t 



FIG. 2: (Color online) Unbiased superdifhision in a wash- 
board potential. The mean square displacement of a particle 
that spreads ballistically in the absence of a bias F, see Fig. 1, 
eventually changes its behavior from superdiffusive to normal 
diffusive behavior under the influence of a periodic poten- 
tial of strength Vo- From bottom right hand side to top we 
use Vo = 0.5, 0.2, 0.1, 0.05, 0.02, 0. The time of the turn-over 
shifts to later times with decreasing potential strength. For 
the simulation of the displayed data the first embedding was 
used with k, = 2, v = 3, ljq = 1 and u(0) = 0. 



C. Superdiffusion in a biased washboard potential 




t 



FIG. 3: (Color online) Biased superdiffusion. After diffusion 
has become normal in the presence of a periodic potential, 
see Fig. 2, it is again changing to ballistic diffusion under 
the influence of an additional finite bias F. Long transients 
exhibiting hyper-diffusion emerge before the ballistic diffusion 
regime is approached. A fixed barrier height Vo = 1 was 
used for the simulations and the tilt F is variable. The other 
parameters are again k = 2, v = 3, ujq = 1 and v(0) = 0. 

The modification of the dynamics by a tilt of the wash- 
board potential, V(x,t) = — Vq cos(27rx/a;o) — Fx, pro- 



vides an intriguing question. In particular one may ask 
whether the spreading will again become superdiffusive 
and whether a supercurrent will emerge that steadily 
grows with time? The numerical simulations displayed in 
Figs. 3 and 4 indicate that the answer to both questions 
is yes. Both the mean square displacement as well as the 
average displacement become proportional to a ballistic 
law t 2 . The time at which this presumably asymptotic 
behavior sets in becomes increasingly larger the smaller 
the bias force F is. For the stronger forces F = 0.7, 1.5 
the ballistic regime has settled within the total time of 
t = 10 4 which requires a week of computational time on a 
Pentium PC with 3 GHz tact-frequency. For small forces 
F = 0.1,0.2 an approach to a ballistic behavior is not yet 
visible. We expect it to occur at a later time. 

Another interesting feature is the occurrence of very 
long superdiffusive transient episodes with a mean square 
displacement growing faster than ballistic as t a " !i with an 
exponent a e g > 2 up to approximately 5. We call these 
episodes "hyper-diffusive" . Their occurrence depends on 
the dimensionless barrier height Vn/fcsT and the biasing 
force i*o- With a larger barrier the total transient time 
before the asymptotic ballistic behavior sets in becomes 
larger. For small biasing forces after a short initial period 
first a regime of normal diffusion is observed which turns 
over into the hyper-diffusive regime at a time that is the 
later the smaller the biasing force is. For example, for 
Vq/UbT — 1 and F = 0.2 the normal diffusion regime 
extends approximately over one decade from t = 10 to 
t = 100, and then rapidly turns around t = 2 x 10 2 into 
hyper-diffusion with a c ff ~ 5.1, cf. Fig. 3. This behavior 
continues until the end of the simulation at t — 10 4 . Until 
then the root mean square displacement increases by an 
amount of 10 3 — 10 4 periods of length x$. The turnover 
to the expected ballistic diffusion can only be observed 
if the biasing force is larger, but then also the normal 
diffusion regime disappears. 

A similar effect of hyper-diffusive motion was reported 
by Lii and Bao [40] for a Brownian particle moving in a 
biased periodic potential under the influence of a super- 
Ohmic model with a spectral density J(w) cx u) for u> — > 
0. The question whether the hyper-diffusion observed 
in Ref. [40] is indeed asymptotic or whether it is also a 
transient phenomenon must still be clarified. 

From the different curves displayed in Fig. 3 one can 
infer that for large times the mean square displacement 
grows the faster the smaller the biasing forces is, in other 
words, the ballistic diffusion constant increases with de- 
creasing biasing force and, in particular, is larger than 
the ballistic diffusion constant of free motion reached for 
Fxo Vo. This phenomenon is akin to the effect of 
giant enhancement of normal diffusion in periodic poten- 
tials [42, 43]. 

The mean displacement (Ax(t)) exhibits a qualita- 
tively similar behavior as the mean square displacement. 
After a first transient period whose nature strongly de- 
pends on the initial velocity distribution a monotonous 
growth sets in that changes from linear to quadratic, pos- 
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sibly interrupted by an episode of rapid growth propor- 
tional to t" with j3 > 2. cf. Fig. 4. The transitions 
between the different regimes occur at the same times at 
which the mean square displacement changes from nor- 
mal diffusion into the hyper-diffusion and finally to bal- 
listic diffusion. The exponent /3 though is much smaller 
than the hyper-diffusive exponent a e s. This indicates 
that the transport in this intermediate regime is strongly 
erratic. While both periods of normal and ballistic diffu- 
sion can be characterized by a time- independent Peclet 
number Pe = x (Ax(t)} / (Ax 2 (t)) [41] the difference of 
the exponents a e ff and /3 does not allow the definition 
of a Peclet number in the hyper-diffusive regime. How- 
ever, both in the normal and the asymptotic ballistic 
regime a time-independent Peclet number can be defined. 
For F = 0.1 the FDT (31) holds with a good accuracy 
and Pe ss Fxq/ (2ksT) in the normal diffusion transport 
regime, see Figs. 3, 4. Beyond the linear response regime, 
the FDT (31) is generally violated. Such a wealth of 
different transport regimes with normal and anomalous 
features, revealed by a simple model is really surprising. 
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clear indication of non-ergodicity. In the case of ballistic 
diffusion in a tilted periodic potential analytic results for 
the velocity fluctuation autocorrelations are not available 
and we therefore have to rely on our numerical findings. 
In Fig. 5 the mean square deviations of position for dif- 
ferent distributions of the initial velocities are compared 
with each other. Apart from minor deviations, these ini- 
tial preparations do not seem to have any influence in 
the presence of a tilted periodic potential. Therefore, 
one might suppose that in this case the process of the 
velocity fluctuations is wide sense ergodic. This result 
though cannot be considered conclusive because also in 
the absence of any potential the choice of the initial dis- 
tribution of velocities has only little impact on the mean 
square displacement, see lines labeled by Vq = 0, F = 
in Fig. 5. A more convincing argument results from the 
comparison of the effect of different initial distributions 
of the auxiliary variables u\ and 112, see Fig. 6. While the 
influence of this distribution on the position mean square 
deviation is very large and even increases with growing 
time, see Fig. 1, only small deviations at early and inter- 
mediate times are visible in the case of a tilted periodic 
potential. Hence, numerical evidence seems to indicate 
that the velocity fluctuations of a ballistically diffusing 
particle in a tilted washboard potential indeed is wide 
sense ergodic. 

This raises the question whether it is possible that the 
velocity fluctuations of a superdiffusive process may be 
wide sense ergodic in one case and non-ergodic in an- 
other. The strict answer to this question is that the 
velocity fluctuations of any truly ballistic diffusion with 
(Ax 2 (t)) — D2t 2 constitute a non-ergodic process. This 
follows from Eq. (23) by means of differentiation with 
respect to time yielding 



FIG. 4: (Color online) Anomalous drift behavior. A finite bias 
F induces anomalous drift. From bottom right hand side to 
top we use F = 0.1,0.2,0.4,1.5. Ballistic currents appear 
asymptotically (Ax(t)} ~ t 2 . Like in Fig. 3 transient regimes 
appear with enhanced particle transport stronger than ballis- 
tic. The used parameters are: Vo = 1, K = 2, v = 3, cjo = 1 
and v(0) = 0. 



IV. ERGODICITY 

We next comment on the ergodic properties of the ve- 
locity fluctuations in relation to ballistic diffusion. In the 
case of free ballistic diffusion, the velocity fluctuations are 
clearly non-ergodic. As already mentioned above, this 
rigorously follows from the fact that the velocity fluctu- 
ation correlation coefficient K v (t) given by Eq. (33) con- 
verges to a constant value different from zero. Also the 
strong dependence of the position mean square displace- 
ment on the initial distribution of the auxiliary variables 
u\ and U2 in the limit of large times, see Fig. 1, provides a 



1 r* 

D 2 = lim - / dt'K v (t') . (36) 

t->oo t J Q 

Therefore the time average of the autocorrelation func- 
tion of the velocity fluctuations does not vanish and con- 
sequently the velocity fluctuations are non-ergodic, see 
Appendix B. 

However, one must keep in mind that the ballistic dif- 
fusion presents a marginal case. Any increase of (Ax 2 (t)} 
slower than t 2 , such as t 2_e with any small, positive e, 
or t 2 1 In t will lead to a vanishing time average of the ve- 
locity fluctuation autocorrelation function by the same 
argument as above. From the numerical point of view 
there is always a limitation of how accurate the scal- 
ing exponent of (Ax 2 (t)} can be determined. Logarith- 
mic corrections are almost impossible to identify. We 
therefore suppose that the observed ballistic diffusion 
in a tilted washboard potential might strictly speaking 
be marginally sub-ballistic and the velocity fluctuations 
wide sense ergodic. 
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FIG. 5: (Color online) Mean square displacement of position 
for different initially distributed velocities. The solid lines 
mark thermally distributed initial velocities (v 2 (0)) = v%, 
whereas dashed lines mark initially zero velocity v(0) = 0. 
The differences in time evolution vanish in the asymptotic 
long time limit in the presence of a biased periodic wash- 
board potential with potential height Vo and bias F, imply- 
ing wide sense ergodicity for the mean square displacement. 
However, in case of free ballistic diffusion, see lines labeled by 
Vo — 0, F — 0, a constant deviation remains according to Eq. 
(35). Parameters are chosen as k = 2, v = 3 and ojq — 1. 



10 8 


F = 1.5 




10 4 






10° 




0.2 


10" 4 







10" 1 10° 10 1 10 2 10 3 10 4 



/ 

FIG. 6: (Color online) Role of deviation from stationary 
fluctuation-dissipation relation in Eqs. (2, A9) on the time 
evolution of the mean square displacement (Ax 2 (t)}. In con- 
trast to the choice with a stationary fluctuation-dissipation re- 
lation, see solid lines marking (ttj(O)itj(O)) = K 2 Sij, the initial 
choice iii(0) =112(0) =0, see dashed lines, yields a Gaussian 
noise f(t) that initially is non-stationary, see Eq. (A8). The 
noise £(t) assumes, however, stationary noise at asymptotic 
long times. The initial velocity was set to zero, i.e. v(0) = 0, 
the potential strength to Vo — 1 and the remaining parame- 
ters are chosen as in Fig. 5, n = 2, v = 3 and u)o = 1- 



V. SUMMARY 

In this work we considered one of the simplest models 
for the superdiffusive motion of a particle described by a 
GLE. It corresponds to a bi-exponential memory kernel 
with zero integral. The according spectral density J(uj) 



of thermal bath oscillators sets in with a cubic law. It 
describes, for example the diffusion of an impurity in a 
crystal. 

We considered a large family of Markovian embed- 
ding schemes, i.e. higher dimensional Markovian pro- 
cesses that generate the considered non-Markovian pro- 
cess upon projection onto the subspace spanned by po- 
sition and momentum of the particle. Out of the whole 
class we identified a simple four-dimensional embedding 
that can be numerically treated in an efficient way. 

We confirm that ballistic superdiffusion is non-ergodic, 
which is concordant with the findings in [25, 26]. As a 
new amazing manifestation of non-ergodicity we found 
that a non-equilibrium initial noise preparation can 
change the law of diffusion, see in Fig. 1. 

Further on our numerical findings indicate that the free 
ballistic diffusion, being present in the absence of any po- 
tential, changes into normal diffusion in the presence of 
a periodic potential. We concluded that the process of 
the velocity fluctuations is non-ergodic in the absence of 
a periodic potential but wide sense ergodic in the pres- 
ence of a periodic potential. Apparently, the transition 
to the ergodic motion does not require a minimal poten- 
tial strength. Rather, the time to reach the asymptotic 
regime of normal diffusion diverges with vanishing poten- 
tial strength Vq. 

An additional biasing force leads to ballistic motion in 
a periodic potential, i.e. both the mean value and the 
variance of the position displacement grow proportional 
to a t 2 law. In this case, however we found strong indica- 
tions that the velocity fluctuations remain wide sense er- 
godic. This paradoxically looking scenario - non-ergodic 
for free ballistic diffusion versus ergodic for ballistic diffu- 
sion in a potential - is possible because the ballistic diffu- 
sion presents a marginal situation. Although for ballistic 
diffusion following a strict t 2 law the velocity fluctua- 
tions are non-ergodic any modification of the t 2 law with 
a weakly decaying function such as 1/ In t leads to wide 
sense ergodic velocity fluctuations. For a sub-critical bias 
F < F = 2nVo I x the ballistic diffusion coefficient D 2 is 
substantially enhanced compared to the diffusion coeffi- 
cient for free ballistic diffusion. This effect is the analog 
to giant enhancement of normal diffusion in tilted wash- 
board potentials. 

Depending on the potential and the bias strengths the 
time before the asymptotic ballistic motion sets in may be 
extremely large. Within this long transient period a nor- 
mal and even a hyper-diffusion regime may exist, where 
a e g exceeds the ballistic value of 2. The presence of such 
long transients presents a general feature of the studied 
non-Markovian dynamics. 

Acknowledgements This work was supported by the 
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Appendix A: Conditions for embedding 

The solution of the last equation in Eq. (7) is 



i(t) = - f 

Jo 



-Mt-t') p ( t r )tdt , + 



/ e- A ^-^C^(t')dt' + e~ At u(0) , (Al) 
Jo 



which inserted into (7) yields 
d 



Pit) 



_dx 
d_ 

dx 



V(x,t)-^(t) 
t 



Te- A{t - t ' ) rp(t')dt'+ 



■V(x,t)-J 

f ^ T e- A ^ci{t')dt' + -$ T e- At H(0) . (A2) 
Jo 

The comparison of (A2) with (6) gives Eq. (9) and 

C(t) = / t ^ T e- A ( t -*')c|(t')^' + 9 T e" At «(0) . (A3) 
Jo 

Assuming (iij(0)£j(i)) = 0, this enables one to calculate 
the noise correlation function: 

<C(*)C(*)> = (9 T e- At ^m T e- A ^(0)) 

+ ( f dt' / S ^ T e- A ( t -*')c|(t')5 T e- A(s - s ' ) Ce(s')) 
Jo Jo 

(A4) 

Taking into account Eq. (8) for t > s (the case t < s can 
be treated alike) the second term reduces to: 

/Ve- A(t - s,) CC T e - AT ( s - s ')^ S ' = 
Jo 

= 5 T e~ A * f S e As 'CC T e ATs 'd S 'e- AT 
Jo 

and the first term is: 

T e- A ^(O)^e^H(O} = 
=^ T e- At @(0)®H T (0))e- ATs ^ 

Altogether, with the following definition: 

U= @{0)®^ T (0} 

this yields: 



(A5) 



(A6) 
(A7) 



— *TV — At 



Ut)((s)) = 9 L e 



U 



s As 'CC T e ATs 'd.s' 



Making an Ansatz as in Eq.(ll) enables one to separate 
the noise correlation function into a stationary and a non- 
stationary part: 



ma*) =Tz 



U + [' e As ' (AG + GA T )e ATs 'ds' 
Jo 



-A r s-> 

x e g 



=-g T e~ At 



U + e As 'Ge ATs ' 



s —s 
s'=0 



-A 1 

e s g = 



^ T e~ A * [U - G] e- AT ^ + Te-^-^G^ 

(A8) 



The first term of Eq. (A8) represents the non-stationary 
part and is vanishing asymptotically in the limit of long 
times, i.e. t, s — > oo. Both the relaxation spectrum defin- 
ing the corresponding time scales and the spectrum of 
autocorrelation times is given by the eigenvalues of ma- 
trix A. Moreover, the fluctuation-dissipation relation of 
Eq. (2) is always asymptotically fulfilled, if one chooses 
G]j = raksT^r , which yields Eq. (10). However, in order 
to obey the fluctuation-dissipation relation for all times, 
one has to set U = G, which implies Eq. (12) and we 
end with: 



(C(t)C(s)) = mk B Tl) T e- A ^-^r = mk B T 1 {t - s) 



(A9) 



Appendix B: Wide sense ergodicity 

According to its definition, a stationary process y(t) 
is ergodic in the wide sense if its time average converges 
in the mean square sense towards the ensemble average. 
This definition implies that a process y(t) is wide sense 
ergodic if and only if the time average of the autocorre- 
lation function of its fluctuations vanishes, i.e. if 



lim — 



dt'((y(t)-(y))(y(0)-(y)))=0 (Bl) 



_ A s— > 

e 9 ■ 



holds [44] . Hence, the decay of the autocorrelation func- 
tion of the fluctuations towards zero provides a sufficient 
condition for a wide sense ergodic process. On the other 
hand, the considered process is non-ergodic if the auto- 
correlation function of its fluctuations approaches a con- 
stant value different from zero. 
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